Spatial Data in R: Mini Project 1

The goals for this phase of your final project are:

Articulate an interesting research question based on a dataset you’d like to learn more about.

Develop a spatial database that contains potentially relevant explanatory variables that you’d like to explore in the context of that research question.

Demonstrate an understanding of the various workflow elements involved in designing and constructing a spatial database for subsequent visualization and analysis.

My Question: Where is the best place to go in New Zealand during the zombie apocalypse to avoid infection and survive? Specifically, I will look at 1. Where the most zombie attacks are happening 2. Are they happening in more urban or rural locations? 3. Highlight rural locations as areas of survival (the less people around, the less chance of encountering an infected person and becoming a zomibe) 4. Number of pharmacies in each city>>region to inform survivors of the best places to go where they will have regular access to medical supplies

I will answer this question using the following data:

library(here)
## here() starts at /home/michaela/R/Mini-Project-1-michaelagustafson
library(tidyr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(terra)
## terra version 1.4.22
## 
## Attaching package: 'terra'
## The following object is masked from 'package:dplyr':
## 
##     src
## The following object is masked from 'package:tidyr':
## 
##     extract
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(tidyverse)
## Warning in system("timedatectl", intern = TRUE): running command 'timedatectl'
## had status 1
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.4     v stringr 1.4.0
## v readr   2.0.1     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x terra::extract()  masks tidyr::extract()
## x dplyr::filter()   masks stats::filter()
## x dplyr::lag()      masks stats::lag()
## x purrr::simplify() masks terra::simplify()
## x terra::src()      masks dplyr::src()
library(nngeo)
library(data.table)
## 
## Attaching package: 'data.table'
## The following object is masked from 'package:purrr':
## 
##     transpose
## The following objects are masked from 'package:terra':
## 
##     copy, shift
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
library(tmap)
## Registered S3 methods overwritten by 'stars':
##   method             from
##   st_bbox.SpatRaster sf  
##   st_crs.SpatRaster  sf
library(viridis)
## Loading required package: viridisLite
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /usr/share/proj
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## 
## Attaching package: 'rgdal'
## The following object is masked from 'package:terra':
## 
##     project
# library(ggmap)

A. Dependent Variable with Point Geometry: Zombie Attacks

Data source: kaggle.com

zom <- read.csv("zombies.csv")

head(zom); dim(zom)
##   zombieid zombie age    sex rurality household water food medication    tools
## 1        1  Human  18 Female    Rural         1     0 Food Medication No tools
## 2        2  Human  18   Male    Rural         3    24 Food Medication    tools
## 3        3  Human  18   Male    Rural         4    16 Food Medication No tools
## 4        4  Human  19   Male    Rural         1     0 Food Medication    tools
## 5        5  Human  19   Male    Urban         1     0 Food Medication No tools
## 6        6  Human  19 Female    Urban         1     0 Food Medication    tools
##                firstaid sanitation clothing documents
## 1    First aid supplies Sanitation Clothing      <NA>
## 2    First aid supplies Sanitation Clothing      <NA>
## 3    First aid supplies Sanitation Clothing      <NA>
## 4 No first aid supplies Sanitation Clothing      <NA>
## 5    First aid supplies Sanitation     <NA>      <NA>
## 6    First aid supplies Sanitation Clothing      <NA>
## [1] 200  14
zombie.shp <- read_sf(here("New_Zealand_Zombie_Attack_Data.shp"))
plot(zombie.shp) #remember to use st_geometry to avoid plotting every column in the dataframe

B. Regions

1. New Zealand Region Data

Data source: https://datafinder.stats.govt.nz/layer/98765-regional-council-2019-clipped-generalised/

The regions of New Zealand are my ‘regions’ I will use to organize higher order variables.

# Find "Regional Council 2018 Clipped (generalised)"
# select the GeoPackage option in the "Vectors/tables" dropdown
# at https://datafinder.stats.govt.nz/data/ (requires registration)
# Save the result as:

unzip("statsnzregional-council-2019-clipped-generalised-SHP.zip")
nz.full.shp <- st_read("regional-council-2019-clipped-generalised.shp")
## Reading layer `regional-council-2019-clipped-generalised' from data source 
##   `/home/michaela/R/Mini-Project-1-michaelagustafson/regional-council-2019-clipped-generalised.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 17 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 1089970 ymin: 4747987 xmax: 2470042 ymax: 6223156
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
plot(st_geometry(nz.full.shp))

# print(object.size(nz_full), units = "Kb") # 14407.2 Kb

# nz <- st_simplify(nz_full, dTolerance = 1000)


names(nz.full.shp)
## [1] "REGC2019_V" "REGC2019_1" "LAND_AREA_" "AREA_SQ_KM" "Shape_Leng"
## [6] "geometry"
nz.full.shp$REGC2019_1
##  [1] "Northland Region"         "Auckland Region"         
##  [3] "Waikato Region"           "Bay of Plenty Region"    
##  [5] "Gisborne Region"          "Hawke's Bay Region"      
##  [7] "Taranaki Region"          "Manawatu-Wanganui Region"
##  [9] "Wellington Region"        "West Coast Region"       
## [11] "Canterbury Region"        "Otago Region"            
## [13] "Southland Region"         "Tasman Region"           
## [15] "Nelson Region"            "Marlborough Region"      
## [17] "Area Outside Region"
nz.reg.shp = filter(nz.full.shp, REGC2019_1 != "Area Outside Region") %>%
        select(Region = REGC2019_1, AreaSQKM = AREA_SQ_KM, geometry)

st_is_valid(nz.reg.shp)
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE
st_crs(nz.reg.shp)
## Coordinate Reference System:
##   User input: NZGD2000 / New Zealand Transverse Mercator 2000 
##   wkt:
## PROJCRS["NZGD2000 / New Zealand Transverse Mercator 2000",
##     BASEGEOGCRS["NZGD2000",
##         DATUM["New Zealand Geodetic Datum 2000",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4167]],
##     CONVERSION["New Zealand Transverse Mercator 2000",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",173,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",1600000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",10000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["New Zealand - onshore"],
##         BBOX[-47.33,166.37,-34.1,178.63]],
##     ID["EPSG",2193]]

C. 1. Census Data

Data source: https://www.stats.govt.nz/topics/census

I will use census data to determine the population density/sqkm in areas of the regions (urban vs rural) and overall population density/sqkm of a region. This information will be best used to determine good regions for individuals to go to be in less populated areas and have a better chance of avoiding other people and potential zombie infection.

# read in census data
nz.census <- read.csv("Census_Usually_resident_population_count_and_change_by_region_2006_2013_and_2018.csv")

# filter to get most recent census population (2018)
names(nz.census)
## [1] "X...Census.year" "Region"          "Region.Code"     "Measure"        
## [5] "Value"           "Value.Unit"      "Value.Label"     "Null.Reason"
nz.census <- filter(nz.census, X...Census.year == "2018") %>%
  select(Year = X...Census.year, Region, Population = Value)

nz.census$Region = paste(nz.census$Region, "Region", sep = " ")
nz.census[8, 2] = "Manawatu-Wanganui Region"

You’ll want to leave yourself a few more “breadcrumbs”, I think. Anytime you are directly modifying the data (as you do above), it’s good to leave a note to remind you why you did

2. Urban vs Rural shapefile

Data source: https://datafinder.stats.govt.nz/layer/98743-urban-rural-2019-clipped-generalised/

unzip("statsnzurban-rural-2019-clipped-generalised-SHP.zip")
urbrur <- read_sf("urban-rural-2019-clipped-generalised.shp")
plot(st_geometry(urbrur))

urban <- urbrur %>%
  select(Name = UR2019_V_1, Urban = IUR2019__1, AreaSQKM = AREA_SQ_KM, geometry)

unique(urban$Urban)
## [1] "Inland water"      "Large urban area"  "Major urban area" 
## [4] "Medium urban area" "Oceanic"           "Rural other"      
## [7] "Rural settlement"  "Small urban area"
# neet to change columsn to say just urban or rural

urban$Urban[grep("urban", urban$Urban)] <- "Urban"
urban$Urban[grep("Rural", urban$Urban)] <- "Rural"
urban$Urban[grep("water", urban$Urban)] <- "Other"
urban$Urban[grep("Oceanic", urban$Urban)] <- "Other"

3. Pharmacy Location Data

Data source: https://fyi.org.nz/request/10183-total-number-of-licenced-pharmacies-in-nz-2019

Does this still count as tabular data even though I turned it into spatial data??

# read in pharmacy csv with location data

nz.pharm <- read.csv(here("nz.pharm.csv"))

nz.pharm <- na.omit(nz.pharm)
nz.pharm.sf <- st_as_sf(nz.pharm, coords = c("lon", "lat"), crs = 4326)
st_crs(nz.pharm.sf)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
nz.pharm.proj <- st_transform(nz.pharm.sf, crs = st_crs(nz.reg.shp))
st_crs(nz.reg.shp) == st_crs(nz.pharm.proj)
## [1] TRUE
plot(st_geometry(nz.reg.shp))
plot(st_geometry(nz.pharm.proj), add = TRUE)

nz.pharm.proj$type <- "Pharmacy"
#looks alright

D. Creating Database

Spatial Join

Questions to Answer: 1. What is the percentage of urban and rural areas within each region? 2. What is the number of attacks in each urban type within each region? 3. What is the population density of each region (people/sqkm)? 4. What is the number of attacks per population in each region? (attacks/people)

# make sure CRS are the same
st_crs(zombie.shp) == st_crs(urban)
## [1] FALSE
zombie.shp.proj <- st_transform(zombie.shp, crs = st_crs(urban))
st_crs(zombie.shp.proj) == st_crs(urban)
## [1] TRUE
st_crs(urban) == st_crs(nz.reg.shp)
## [1] TRUE
# ended up using a nearest neighbor approach because I kept getting points 800 something rows when I started with 668 so something was getting duplicated:
#zombie.urban.join <- st_join(urban, zombie.shp.proj, join = st_nn, k = 1, maxdist = 50)


#zombie.urban.region.join <- st_join(zombie.urban.join, nz.reg.shp, join = st_within)


#str(zombie.urban.region.join)

#zombie.urban.region.ng <- st_drop_geometry(zombie.urban.region.join)

#zombie.urban.region.df <- zombie.urban.region.ng %>%
 # select(Urban, AreaSQKM.x, Attacks, Region)



#zombie.urban.region.df <- as.data.frame(zombie.urban.region.df)

#zombie.urban.region.df$Region <- as.factor(zombie.urban.region.df$Region)

#zombie.urban.region.df$Attacks <- replace_na(zombie.urban.region.df$Attacks, 0)
#zombie.urban.region.df$Attacks <- as.numeric(zombie.urban.region.df$Attacks)


#summary.df <- zombie.urban.region.df %>% 
  #group_by(Region, Urban) %>% 
  #summarise(., sum(Attacks), sqkm_dis_total = sum(AreaSQKM.x))

#summary.df <- rename(summary.df, Attacks = "sum(Attacks)")

#summary.df$urban_attacks <- paste(summary.df$Urban, summary.df$Attacks, sep = ";")

#summary.df.wide <- spread(summary.df, Urban, sqkm_dis_total, "0")
#summary.df.wide <- spread(summary.df.wide, urban_attacks, Attacks)

#str(summary.df.wide)

#summary.df.wide$Other <- as.numeric(summary.df.wide$Other)
#summary.df.wide$Rural <- as.numeric(summary.df.wide$Rural)
#summary.df.wide$Urban <- as.numeric(summary.df.wide$Urban)

#summary.df.wide <- summary.df.wide[!is.na(summary.df.wide$Region),]

#summary.df.wide[is.na(summary.df.wide)] <- 0

#summary.df.wide <- summary.df.wide %>%
 # group_by(Region) %>%
  #summarise_if(., is.numeric, sum)


#summary.df.wide <- rename(summary.df.wide, Other_sqkm = "Other")
#summary.df.wide <- rename(summary.df.wide, Rural_sqkm = "Rural")
#summary.df.wide <- rename(summary.df.wide, Urban_sqkm = "Urban")
#

#colnames(summary.df.wide)[grepl('Other;',colnames(summary.df.wide))] <- 'Other_attacks'
#colnames(summary.df.wide)[grepl('Rural;',colnames(summary.df.wide))] <- 'Rural_attacks'
#colnames(summary.df.wide)[grepl('Urban;',colnames(summary.df.wide))] <- 'Urban_attacks'

#names(summary.df.wide) <- make.unique(names(summary.df.wide))
#tidy.summary.df <- summary.df.wide %>%
  #rowwise() %>%
  #mutate(Rural_Attack_sum = sum(across(starts_with("Rural_attacks")), na.rm = T))

#tidy.summary.df <- tidy.summary.df %>%
  #rowwise() %>%
 # mutate(Other_Attack_sum = sum(across(starts_with("Other_attacks")), na.rm = T))

#tidy.summary.df <- tidy.summary.df %>%
  #rowwise() %>%
  #mutate(Urban_Attack_sum = sum(across(starts_with("Urban_attacks")), na.rm = T))

# taking out Other Attack sum becuase there were none (they were in an unidentified region)
#tidy.summary.df <- tidy.summary.df %>%
  #select(Region, Other_sqkm, Rural_sqkm, Urban_sqkm, Rural_Attack_sum, Urban_Attack_sum)
### NEED TO FIX REGIONAL SQKM
# Find the number of pharmacies per region

pharm.per.region <- st_join(nz.reg.shp, nz.pharm.proj, join = st_contains)
pharm.per.region.ng <- st_drop_geometry(pharm.per.region)
pharm.per.region.ng <- pharm.per.region.ng %>%
  select(Region, Pharmacy_Nm)

# Count the number of pharmacies in each region
pharm.per.region.count <- pharm.per.region.ng %>% 
  group_by(Region) %>%
  tally()

pharm.per.region.count <- rename(pharm.per.region.count, pharms_per_region = n)
# combine pharmacies per region with summary.df

#tidy.summary.df <- left_join(tidy.summary.df, pharm.per.region.count)

Add census data to tidy summary df

# add total population to final df
#tidy.summary.df <- left_join(tidy.summary.df, nz.census)

# add total area (sqkm) to final df
#tidy.summary.df <- left_join(tidy.summary.df, nz.reg.shp)

# calculate attacks/rural area and attacks per urban area per region

#tidy.summary.df <- tidy.summary.df %>%
  #mutate(attacks_per_rural = Rural_Attack_sum/Rural_sqkm,
        # attacks_per_urban = Urban_Attack_sum/Urban_sqkm,
        # pop_density = Population/AreaSQKM,
        # perc_rural = Rural_sqkm/AreaSQKM,
        # perc_urb = Urban_sqkm/AreaSQKM,
        # total_attacks = Rural_Attack_sum + Urban_Attack_sum)

<<<<<<< HEAD # Add elevation

rastlist1 <- list.files(path = "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/", pattern = '.tif', full.names=TRUE)
rastlist1
## [1] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//aklnd_25r.tif"   
## [2] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//ecape_25r.tif"   
## [3] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//hksbay_25r.tif"  
## [4] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//nthcape_25r.tif" 
## [5] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//taranaki_25r.tif"
## [6] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//waikato_25r.tif" 
## [7] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//well_25r.tif"    
## [8] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//whngrei_25r.tif"
rn1 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/aklnd_25r.tif")
rn2 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/ecape_25r.tif")
rn3 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/hksbay_25r.tif")
rn4 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/nthcape_25r.tif")
rn5 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/taranaki_25r.tif")
rn6 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/waikato_25r.tif")
rn7 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/well_25r.tif")
rn8 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/whngrei_25r.tif")


rastlist2 <- list.files(path = "/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/", pattern = '.tif', all.files=TRUE, full.names=FALSE)
rastlist2
## [1] "chch25r.tif"     "dunedin_25r.tif" "greymth_25r.tif" "invcgll_25r.tif"
## [5] "kaik_25r.tif"    "mtcook_25r.tif"  "nelson_25r.tif"  "teanau_25r.tif" 
## [9] "waitaki_25r.tif"
rs1 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/chch25r.tif")
rs2 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/dunedin_25r.tif")
rs3 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/greymth_25r.tif")
rs4 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/invcgll_25r.tif")
rs5 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/kaik_25r.tif")
rs6 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/mtcook_25r.tif")
rs7 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/nelson_25r.tif")
rs8 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/teanau_25r.tif")
rs9 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/waitaki_25r.tif")

all.rast <- merge(rn1, rn2, rn3, rn4, rn5, rn6, rn7, rn8, rs1, rs2, rs3, rs4, rs5, rs6, rs7, rs8, rs9)
 
zombie.elev.proj <- st_transform(zombie.shp, crs = crs(all.rast))

plot(all.rast)
plot(zombie.elev.proj, add = TRUE)
## Warning in plot.sf(zombie.elev.proj, add = TRUE): ignoring all but the first
## attribute

# Technically should be calculating the mean elevation or range or something from a buffer surrounding the point, but for the sake of time (because I already tried it and I'd have to mess around with changing extents and a bunch of other stuff) and to actually practice what I'm supposed to be practicing - just going to the pull the elevation from the exact point of attack:

zomb.elev.proj.vect <- vect(zombie.elev.proj)
elev.at.attack <- extract(all.rast, zomb.elev.proj.vect, na.rm = TRUE)

elev.at.attack$ID <- zombie.shp$OBJECTID
elev.at.attack <- rename(elev.at.attack, OBJECTID = ID)

Trying to create new df

# add urban or rural  
zombie.attacks.df <- st_join(zombie.shp.proj, urban, join = st_within)
colnames(zombie.attacks.df)
##  [1] "OBJECTID"   "OBJECTID_1" "Day"        "Attacks"    "Time"      
##  [6] "Longitude"  "Latitude"   "geometry"   "Name"       "Urban"     
## [11] "AreaSQKM"
zombie.attacks.df <- zombie.attacks.df[,c(1, 4, 6:8, 10)]

# add elevation
zombie.attacks.df <- left_join(zombie.attacks.df, elev.at.attack)
## Joining, by = "OBJECTID"
zombie.attacks.df <- rename(zombie.attacks.df, elev = aklnd_25r)

# add region
zombie.attacks.df <- st_join(zombie.attacks.df, nz.reg.shp, join = st_within)

# add number of pharmacies per region
zombie.attacks.df <- left_join(zombie.attacks.df, pharm.per.region.count)
## Joining, by = "Region"
# add census data
zombie.attacks.df <- left_join(zombie.attacks.df, nz.census)
## Joining, by = "Region"
# mutate and calculate population density per region
zombie.attacks.df <- zombie.attacks.df %>%
  mutate(pop_density = Population/AreaSQKM)

#####

nz.reg.shp.join <- left_join(nz.reg.shp, nz.census)
## Joining, by = "Region"
nz.reg.shp.join <- nz.reg.shp.join %>%
  mutate(pop_density = Population/AreaSQKM)

======= you didn’t have any raster extractions here and the analysis is not reproducible because a) you didn’t save the nz_pharm object, b) you didn’t provide all of the necessary info to recreate it. I would have liked to see a little more effort and annotation. If I knew more about what you were trying to do, I might actually be able to fix the parts of this that wont run.

save.image(here("workspaces/MiniProj_Workspace.RData"))

Mini Project 2

library(tidyverse)
library(pander)
## 
## Attaching package: 'pander'
## The following object is masked from 'package:terra':
## 
##     wrap
library(sf)
library(units)
## udunits database from /usr/share/xml/udunits/udunits2.xml
library(ggmap)
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
## 
## Attaching package: 'ggmap'
## The following object is masked from 'package:terra':
## 
##     inset
library(cartogram)
## 
## Attaching package: 'cartogram'
## The following object is masked from 'package:terra':
## 
##     cartogram
library(patchwork)
## 
## Attaching package: 'patchwork'
## The following object is masked from 'package:terra':
## 
##     area
library(tmap)
library(viridis)
library(tigris)
## To enable 
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.
library(ggspatial)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggmap':
## 
##     wind
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

Static Map 1: Location of attacks and population density

zombie.attacks.df <- st_as_sf(zombie.attacks.df)
?tm_polygons

tm_shape(nz.reg.shp.join) +
  tm_polygons(col = "Region", palette = "viridis") +
  tm_legend(outside = TRUE) + 
  tm_bubbles(size = "Population", col = "pop_density", border.col = "black")
## Some legend labels were too wide. These labels have been resized to 0.66. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## The legend is too narrow to place all symbol sizes.

# required libraries
library(leaflet, quietly = T, warn.conflicts = F)

# start basemap
map <- leaflet() %>% 
  
  # add ocean basemap
  addProviderTiles(providers$Esri.OceanBasemap) %>%
  
  # add another layer with place names
  addProviderTiles(providers$Hydda.RoadsAndLabels, group = 'Place names') %>%
  
  # add graticules from a NOAA webserver
  addWMSTiles(
    "https://gis.ngdc.noaa.gov/arcgis/services/graticule/MapServer/WMSServer/",
    layers = c("1-degree grid", "5-degree grid"),
    options = WMSTileOptions(format = "image/png8", transparent = TRUE),
    attribution = NULL,group = 'Graticules') %>%
  
  # focus map in a certain area / zoom level
  setView(lng = 174, lat = -41, zoom = 5)%>%
  
  # add layers control
  addLayersControl(overlayGroups = c('Place names',
                      'Graticules',
                      'Points',
                      'Lines',
                      'Polygons'),
                      options = layersControlOptions(collapsed = FALSE),
                      position = 'topright') %>%
      
      # list groups to hide on startup
      hideGroup(c('Place names'))

# show map
map
### Add data

nz.shp.reproj <- st_transform(nz.reg.shp.join, crs = st_crs(zombie.shp))

library(htmltools)
## 
## Attaching package: 'htmltools'
## The following object is masked from 'package:pander':
## 
##     p
labs <- as.list(zombie.shp$Attacks)
labs2 <- as.list(nz.shp.reproj$pop_density)
## Color setup
colpal <- colorNumeric(palette = "magma", domain=nz.shp.reproj[['pop_density']], n=10)
colorData <- nz.shp.reproj[["pop_density"]]

map <- map %>%
  addTiles() %>%
    addPolygons(data = nz.shp.reproj, label = lapply(labs2, HTML), color = ~colpal(colorData))

map <- map %>%
  addCircleMarkers(data = zombie.shp,
                   weight = 0.5,
                   label = lapply(labs, HTML),
                   col = 'black', 
                   fillColor = 'darkslategrey',
                   radius = 4, 
                   fillOpacity = 0.9, 
                   stroke = T)



map

Honestly, just pretty stoked I got an interactive map to work. When you hover over points and shows how many attacks occured at that point and when you hover over regions it shows the population density (people/sqkm)

These Are My Tears